Existence and uniqueness results for the Gradient Vector Flow 
and geodesic active contours mixed model 

Laurence GUILLOT*, Maitine BERGOUNIOUX* 
February 2, 2008 



*UMR 6628-MAPMO, Federation Denis Poisson, Universite d 'Orleans, BP. 6759, F-45067 
Orleans Cedex 2, laurence.guillot@univ-orleans.fr, maitine.bergounioux@univ-orleans.fr 

Abstract 

This article deals with the so called GVF (Gradient Vector Flow) introduced by C. Xu, 
J.L. Prince [14j[T5]. We give existence and uniqueness results for the front propagation 
flow for boundary extraction that was initiated by Paragios, Mellina-Gottardo et Ralmcsh 
[TT| [12] . The model combines the geodesic active contour flow and the GVF to determine 
the geometric flow. The motion equation is considered within a level set formulation to 
result an Hamilton-Jacobi equation. 

Keywords: image segmentation, gradient vector flow, geodesic active contour, Hamilton- 
Jacobi equation, viscosity solution. 
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1 Introduction 

We consider an image segmentation model : I is a given image and we want to detect bound- 
aries without connexity or convexity assumptions on contours. Therefore we are interested 
in the Gradient Vector Flow (GVF) as a front propagation flow model. This model builds 
a class of vector fields derived from images and has been introduced by Chenyang Xu and 
Jerry L. Prince in [14J. The GVF can be viewed as external forces for active contour models: 
it allows to solve problems where classical methods convergence fail to deal with boundary 
concavities. On the other hand, a new front boundary-based geometric flow for boundary ex- 
traction was proposed by Paragios, Mellina-Gottardo and Ralmesh in 112]. In this model 
the GVF is used to revise the geodesic active contour model of V. Caselles, R. Kimmel, G. 
Sapiro [2] resulting on a bidirectional geometric flow. The classical parametric active contour 
model was proposed by D.Terzopoulos, A.Witkin et M.Kass [13]. It derives from an energy 
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functional minimization. Curves are drawn toward the object boundaries under potential 
external forces action which can be written as the negative gradient of a scalar potential 
function derived from images. Other forces (as pressure forces) may be added. However, the 
performance of methods is limited by unstable initialization process and poor convergence 
when the boundary is concave. Chenyang Xu and Jerry L. Prince [H] set a new external 
force, dealing with these limitations. 

Most of time snakes model provide a local minimum of the functional cost. So, C.Xu 
et JL. Prince generalized this model from the balance equation at the equilibrium between 
internal and external forces (Euler equation of the minimization problem) . They replace the 
standard external force F^l = ~~ VP X (where P is a potential edges detector), by a more 
general external force. 

This new external force field is called the Gradient Vector Flow (GVF). It includes a 
divergence- free component and a curl-free component 116] . Therefore, this new active 
contour model cannot be formulated as an energy minimization problem. The external force, 
denoted V below is introduced via the balance equation that can be written : 

The resulting parametrized curve solving the equation (jl.ip is called " GVF-snake " . The 
final configuration of a GVF-snake satisfies an equilibrium equation which is not a variational 
problem Euler equation, since V(x, y) is not an irrotational field. 

In next section, we introduce the GVF and present the advantages of such a field in 
view of an active contour model and we present and give a precise definition of the Gradient 
Vector Flow. In section 3. we present this model which is inspired by the geodesic active 
contour model combined with the GVF. Finally, the level set strategy leads to an Hamilton- 
Jacobi equation that has not been studied yet (to our knowledge): in the last section we give 
existence and uniqueness results. 

2 The Gradient Vector Flow (GVF) 

The GVF V is a 2-dimensional vector field that should minimize the following objective 
function [13 Q3]. 

E(V) := / n(ul + u 2 y + v 2 x + v 2 ) + f\Vf\ 2 \V - V/| 2 dx (2.1) 
Jn 

where Q is an open, bounded subset of R 2 , u x ,u y , v x ,v y denotes the spatial derivatives (with 
respect to x and y) of the field V = (u, [i > and / : f2 — > R a continuous edge detector. 
There are many choices for /: in [15] . C.Xu and JL. Prince consider 

f(x,y) = Ei x \(x,y) = -\VI(x,y)\ 2 
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where I is the (intensity) image or 

f(x,y) = E®(x,y) = -\V[G a * I{x,y)}\ 2 

where G a is a Gaussian kernel, so that the image is filtered. Here we choose detector proposed 
by R. Deriche and 0. Faugeras in [3] 

f(x,y) = h(\V(G (T *I)(x,y)\ 2 ) 

where 

h(b) = l--^e-& . (2.2) 
V lira 



Here | • | denotes the M 2 - euclidean norm. 

The first term of the functional E is a regularization term whereas the second term is a 
data-driven component. If /|V/| 2 is small, the energy is dominated by the first term and we 
get a slowly varying field. On the hand, when f\Vf\ 2 is large, the second term forces V to 
decrease to V/. Therefore V is close to the gradient of the edge map when it is large (this 
is classical external force for snakes) and does not evolve quickly in homogeneous regions 
(which increase the snake capture area). The parameter \i govern the tradeoff between the 
two integrands of the cost functional and should be set according to the image noise level 
(more noise increases \x) see |15j . 



2.1 On GVF existence 

First we precise the notations : Q is an open, bounded subset of M. 2 with C°° boundary T: 
it is the image domain. Let V = (u,v) and W = (£, x) m H(fi) : = H 1 ^) x H 1 ^); then 
DV(x) = (V«(x),Vw(x)) G L 2 (Q) 4 . The inner product in L 2 (Vtf is 

(DV,DW) 2 = J (Vt*,VO*» + (Vv,Vx)Ridx = (Vu, V£) i2(n) 2 + (Vv,V x ) L 2 {Q f, 

and the L 2 (ft) 4 -norm of DV is denoted ||W|| 2 . The space W(fi) = H 1 ^) x JT 1 ^) is 
endowed with the norm 

\\Vf n u) = \\V\\ 2 L2{nf + \\DV\\l 
A first "definition" of the GVF could be the following: 

Definition 2.1 The Gradient Vector Flow field is defined as a solution of the following op- 
timization problem 

(V) : min{E(V) | V G (2.3) 

where the energy functional E(V) is defined by h2.1\) . Unfortunately, such a definition is not 
correct since the minimum is not necessarily attained. The study of E will provide another 
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definition of the GVF as the solution to a decoupled system parabolic partial differential 
equations. 

From now and in the sequel we assume that the edge function / verifies (Hi) : 

(Hi) : f G C 1 ^) and / > 

Remark 2.1 Note that if f satisfies (Hi) then f G H 1 ^) and /|V/| 2 G L°°(ft) which are 
the "minimal" assumptions in a first step. 

Proposition 2.1 : Assume (Tti). The functional E is continuous onTL(Vi). 
Proof - We get 

E(V)=»\\DV\\ 2 + [ /|V/| 2 (|Vf -2<y,V/) M2 ) ( fx+ / f\Vf\ 4 dx = TT(V,V) + L(V) + C 
Jn Jn 

tt(X, Y) = fi(DX, DY) 2 + / f\Vf\ 2 (X, Y) R2 dx, (X, Y) G (H(n)) 2 , 

Jn 

L(X) = -2 [ f\Vf\ 2 (X,Vf)dx, X G H(Sl) and C = [ /|V/| 4 dx . 
Jn Jn 

The bilinear form it is continuous on TL(VL) for the i7 1 (ri) 2 -norm : let X et Y be in Tt(Q,), 



where 



\tt(X,Y)\ 



H(DX,DY) 2 + [ f\Vf\ 2 (X,Y)dx 
Jn 

</xpx|| 2 py|| 2 + ||/|v/| 2 ||^ (Q) ||x|| L2(Q)2 ||y|| L2(n)2 
< f)\\ x \\m(si) 2 \\ Y \\m(n) 2 

where C(fx, f) is a constant that only depends on data. L is obviously linear and continuous 
on 7i(VL). We deduce the continuity of E. □ 

Theorem 2.1 The functional E is Gateaux- differentiate on Tt(Q) and for every V = (u, v) 
and W = (£, %) in TL(VL) we get 

(VE(V), W) = 2 f n(DV, DW) + f\Vf\ 2 (V - V/, W)dx (2.4) 
Jn 

Furthermore, if fi G R + then E is convex on 7i(fl). 

Proof - The Gateaux -differentiability is clear. In addition for every V = (u, v) and W = (£, \) 
in H(£l) we have 

E(V + W)- E(V)- (VE(V),W) = [ fi\DW\ 2 + f\Vf\ 2 \W\ 2 dx , 

Jn 
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so that, if n > 0, as / > by (Hi) E(V + W)- E(V) - (VE(V),W) > 0. So we get the 
convexity of E on TC(Q). □ 

From now we assume \i > 0. 

Theorem 2.2 Assume that /|V/| 2 does no£ degenerate on f2 : 

3c> such as f\Vf\ 2 > c (2.5) 

then 7T is coercive on TL{Q). This implies the coercivity and the strict- convexity of E on 
and the problem (V) has a unique solution. 
Proof - it has been define above and we get 



tt(X,X) =fi\\DX\\l+ [ f\Vf\ 2 \X\l 2 dx> n\\DX\\l + c [ \X\l 2 dx > min(/i, c)\\X\\ Hl 

Jn Jn 



(nf 



therefore tt is H 1 ^) -coercive on H(O) and strictly convex. Therefore E is coercive and 
strictly convex on H.(Q) as well. It follows that if (|2.5p is verified then (V) admits a unique 
solution. □ 
Let us formally write the Euler-equations of problem (V) : assume V* = (u* ,v*) is a 
solution to (V). By convexity, it is a stationary point and VE(V*) = 0. Let W = (£, x) S 
H(Q). Integrating by parts expression (|2.4I ) gives 

(X7E(V*),W) 

= -pji&u* + X At;*)dx + J f\Vf\ 2 (V* - V/, W) R2 dx + J + do- (2.6) 

where ^(x) is the outer unit normal of d£l = T at x. 

We first suppose that £ and x belong to V(Q) (the space of C°°(0) functions with compact 
support in Q). Then 



So 



/ (£Au* + X Av*)dx + [ f\Vf\ 2 (V* -Vf,W) R 2dx = 
Jn Jn 



j -fiAu* + f\Vf\ 2 (u* -f x ) = in V(Q) 

1 -fiAv* + f\Vf\ 2 (v* - /„) = in V(Sl) 
where f x , f y are spatial derivatives of /. With (|2.6p . we obtain : 



( 



/ du* \ 






!■( 






\ du J 





el , r I) -u 
xlr 
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Finally, if a solution V* = (u*,v*) to problem (V) exists it must verify : 



-fiAu* + f\Vf\ 2 (u* -f x ) = in n, — = on T 



-vAv* + f\Vf\ 2 (v*-f y ) = inQ, 



Dv 

dv* 



(2.7) 



on r 



The above equations are equilibrium equations if V* realizes the minimum of energy E. 
However, we cannot ensure the existence of such a minimum. Indeed, assumption (|2.5p is not 
realistic : the same gray level for image I on a significant area implies V/ = 0. In this case 
the functional E is a priori non longer coercive and we do not know if (V) has a solution. 
So, instead of computing the "exact" GVF, minimizer of the functional E, we approach it by 
a minimizing sequence and we consider it is the equilibrium state of a time evolving vectors 
fields. The stationary problem becomes a dynamic one : 



V(x) 



lim V(t,x). 



dV 



A simple way to impose a motion to the vectors field is to impose the velocity — — (x, t) 

ot 



setting : 



dV 
~dt 



(x f i) = -V v B(V(i 1 x)) 



(2. 



This leads to parabolic partial differential equations. The gradient vector flow is initialized 
as the gradient of the edge detector / : 

V(t = 0, x) = Vb(x) = V/(x) sur O (2.9) 

We obtain the following dynamic formulation which is the GVF suitable definition : 

Definition 2.2 The gradient vector flow V = (u,v) is defined by the following decoupled 
equations respectively verified by each of its coordinates u and v : 



du 



dv 

{ u(0, •) = / : 



fiAu- f\Vf\ 2 (u- f x ) 




m]0,+oo[x$7 

on ]0, +oo[xr 
in 



(2.10) 



r dv 

ov 

{ v(0,-)=fy 



liAv- f\Vf\ 2 (v - f y ) = in]0,+oo[xft 

on ]0, +oo[xr 
in Vt 



(2.11) 
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2.2 GVF properties. 

In this subsection, we give regularity properties of solutions to (|2.10p and (|2.11|) . Of course, 
it is sufficient to consider equation (|2.f0p . First, the existence of a unique solution is given 
by a classical theorem (see for example [HI HIE])- The bilinear form a associated to equation 
(12301) is 

a(t;u,v) = n (Vu(x).Vv(x) + u(x)v(x)) dx + / f\ V f\ 2 u(x)v(x)dx 
Jn Jo, 

where fx > and /|V/| 2 G L°°(]0, T[xfi).It satisfies : 

1. 1 1— ► a(t;u,v) is measurable Vu, t> £ 

2. For almost i G [0, T] and for all u, v G fl^fi), 

|a(t;u,v)| < C(^, /)||n|| H i (n ) ||w||jfi(n) 

3. For almost t G [0, T] and for all v G JET^fi), 

afcu.t;) > Mll«ll^i(n) + Wf\ V f\ 2 \\L~>(n)\v\ 2 L 2 {n) 

Note that here a does not depend on t. So we may assert that ()2.10j) has a unique solution 
u G L 2 (0,T;H 1 (Q))nC([0,T};L 2 (n)) and — G L 2 (0, T; i? 1 ^)'). Moreover 

Theorem 2.3 Let T > 0, and assume (Hi). Then the GVF (u,v) (solution of $2.1(J\) and 
(TOT]) ) is C 1 on [0,T] x H. 

Proof - We use a generic regularity result ([9]). We prove the result for the component u. 
Assumption (Hi) yields that V/ G L 2 (n) 2 and / |V/ | 2 V/ G [^([O.T] x_H)] 2 . So the solution 
u of ( (12301) belongs to ^([e, T) x H) for all e > 0. Moreover V/ G C°(Q) according to (Hi) 
and compatibility conditions are satisfied with respect to boundary and initial data. So we 
may conclude. □ 
The GVF is built as a spatial diffusion of the edge detector / gradient. This is equivalent 
to a progressive construction of the gradient vector flow starting from the object boundaries 
and moving toward the flat background. In [15], the GVF is normalized to obtain a more 
efficient propagation. It is denoted V^(x) = (u(x),?)(x)) where 

u(x) = =, v(x) 



I 2 2" / 2 ~ 

yit(x) +f(x) Y u ( x ) + ^(x) z 



to give the new external force of the geometric flow called GVF-snake model. The velocity 
of the contour C is given by the equation : 

C t (s,t) = a^-p^ + V(x) (2.12) 
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C. Xu and J.L. Prince have shown that a such flow is not dependent of initial conditions and 
deal with concave regions. However, it depends on curve parametrization, cannot manage 
topological changes, and involves second and fourth order derivatives that are difficult to 
estimate. The objective of N. Paragios, O. Mellina-Gottardo and V. Ralmesh in [llj is to 
eliminate these shortcomings by integrating the GVF with the geodesic active contour and 
implementing it using the level set method [12] . Our aim is to study the Hamilton- Jacobi 
equation derived from this model. 



3 Geodesic contours and GVF 

3.1 Paragios - Mellina-Gottardo - Ralmesh model 

N. Paragios, O. Mellina-Gottardo and V. Ralmesh have defined in [TT] a new "front prop- 
agation flow for boundary extraction". Their geometric model is inspired by the geodesic 
active contours ([2]) and directly defined by the contour evolution velocity. It is based on 
the remark that the Gradient Vector Flow field after the rescaling refers to the direction that 
has to be followed to locally deform the contour and to reach the closest object boundaries. 
On the other hand, given the fact that the propagation of a contour often occurs along the 
normal direction, the propagation will be optimal when V and the unit outward normal v 
are colinear. So we choose to project the normalized gradient vector flow onto the outward 
normal. Then we multiply the velocity by an edge detector function g (that may be different 
from /), which represents the contour information. The contour evolution velocity is then 
given by the equation : 

Ct(x) = ff(|V^(x)| 2 )J(n,i))(x),^(x)) R2 ^ (x) (3.13) 

boundary projection 

where I a is the filtered image and g(b) = ^ e 2^ = 1 — h(b) where h is defined by (|2.2[) 

When there is no boundary information ( |V/<j| 2 <C 1), the contour evolution is driven by 
the inner product between the Normalized Gradient Vector Flow (NGVF) and the normal 
direction : it is adapted to deal with concave regions. When the curve reaches the object 
boundaries neighbourhood ( |V/ CT | 2 ~ +oo ) then g ~ that is, the flow becomes inactive 
and the equilibrium state is reached. 

It is classical to impose a regularity condition on the contour propagation adding a cur- 
vature term and a "balloon force" H. The evolution equation becomes 

C t (x)= 5 (|V/(x)| 2 ) I -^(x)+(l-| J H"(x)|)((u^)(x),i/(x)) R 2+ F(x) |i/(x) (3.14) 



smoothness boundaries attraction balloon force/ 



where (3 > 0. 
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3.2 Level set implementation 



Problems of topologic changes can be solved using the level set method [10]. The moving 2D- 
curve is viewed as the zero level set of a 3D surface which equation is z — = 0. We 

denote $j the partial derivative — — of <!> towards t. 

at 

Theorem 3.1 The evolution of the 3D-surface is described by : 
( * t (x) =< 7 (|V/ ff (x)| 2 )((/3 K (x)- J ff(x))|V$(x)| 



-(l-|^(x)|)(F(x),V$(x)) I 



<9$ 

ov 



\ *(<V) = <*>o 

where <I> is the signed distance defined by : 

{ Oi/xer, 



in ]0, +oo[xS7, 
on ]0, +oo[xT, 
in 



(3.15) 



^o(x) 



± d the distance between x and T 



the positive sign (resp. negative) is chosen if the point x is outside (resp. inside) V. 

Proof - The curve C is the zero set level of <3? : $>(t,C(x)) = 0. Deriving formally with respect 

to t gives for almost every x £ fi: 

$ t (x) + (V$(x),Ct(x)> R2 = 
$t(x) + <V<J?(x), ^(IV/^x)! 2 ) f-^(x) + (1 - |#(x)|)(t>(x),i/(x)) M 2 + (x)) i/(x)) Ra = 



. . V$(x) . , . , . , . . 

where i/(x) = M is the outward unit normal. A short computation gives 



|V*(x) 

* t (x) = 5(|V/ CT (x)| 2 ) ( Wx) - (1 - |i?(x)|)(t>(x), 



Finally 



V$(x) 
|V$(x) 



<I> t (x) = < 7 (|V/ (7 (x)| 2 ) 



#(x) ) |V*(x) 
\ 



(£«W -£(x))|V$(x)| - (1 - |g(x)|) (y(x),y$(x)) 

v^T " (b) 



□ 



The final flow can be decomposed in 

• (a) : a term that provides propagation regularity aims and shrinks the curve toward the 
object boundaries, 

• (b) : a bidirectional flow that moves the curve toward the internal and external objects 
boundaries, 

• (c) : an adaptative balloon force that drives the propagation of the curve when the boundary 
term becomes inactive 1121. 
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3.3 Choice of the edge detector / 

As mentionned before, we focus on the Gaussian edge detector proposed by R. Deriche et O. 
Faugeras in [3]: 

f(x,y) = h(\V(G a *I)(x,y)\ 2 ) (3.16) 

where h is defined by (|2.2p . We must ensure that assumption (Hi) is satisfied for the edge 
detector /. Therefore, we have to set additional hypothesis the image intensity function that 
we have denoted /. 

/ is supposed to have compact support included in f2 that is image frame: we decide 
that I is equal out of the "true image". I G L°°(n) C L 2 (9) C L x (0). On the other 
hand / <£ H 1 ^). Indeed, the image gradient norm becomes infinite at objects contours. 
Furthermore / ^ C 1 (Q) . We correct this lack of regularity using a filtering process that 
makes the filtered image C°° (and of course C 1 ). 

The filtered image is supposed to be C°°, with compact support. So, we cannot choose 
G a * I since the resulting image has no compact support. So we consider a fixed compact 
subset X of Vt and Y a compact subset of fi containing X. 
Let us consider G^ a C°°(Q) projected of the Gaussian kernel G a such as : 

r>*M - I if x e X 

^ CT W ~ \ if x i Y 

So the regularized (filtered) image I a = G* * I verifies 

l a eC?(n) 

and 

Vx G n, V/ ff (x) = V(G* * J)(x) = VG* * /(x). 



Of course the support of I a , contained in supp(G^) + supp(J), is not necessarily included in 
n,, but even if we must extend the frame f2 of the image, we consider that the filtered image 
has a compact support in f2. More precisely : 

Proposition 3.1 The function 



x^V/ CT (x) (3 ' 17) 
is C°°(0,), with compact support in f2 and so bounded on Q. 

Via g L°°(n) nc c °°(0) 

From now and in the sequel, we denote I a by I and we make the following hypothesis on I: 

(Hi) : I eC^ttynW 1 ' 00 ^) 
Note that / may be extended by to £1. Now we can give / regularity properties : 
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Proposition 3.2 Assume the image I satisfies (Hi). The function f defined by : 



x ^ /(x) = M| VJ(x)| 2 ) = 1 " ^^e- — 



belongs to C°°(Q). 

Proof - It is clear since I is C°°(f2), with compact support in f2 and /i is obvioulsy C°°. In 
addition / is nonnegative and bounded by 1 on Q. □ 
Let us denote g the function 

xGSIh 5 (|V/(x)| 2 ) = -^^e S^ - 



27ra 

We have seen in Proposition 13.21 that g = 1 — / G C°°(Q). 
In the sequel, we need the following lemma : 



Lemma 3.1 The function \fg belongs to ; it is Lipschitz continuous on f2 with con- 

stant K\. 



4 Propagation equation study 

Equation (|3. 15|) has been obtained quite formally with the level-set method. Now, we give 
existence and uniqueness of a viscosity solution (p] for example) . Of course, we assume 
that the image I verifies (Hi) so that / satisfies (Hi). The outward unit normal v is a C 1 ' 
Lipschitz vector field. 

4.1 Viscosity theory framework 

Let us precise the notations: forn= ( ^ l ) G M 2 the matrix —pr [ ^ 1 ^ 1 f 2 ] is denoted 

V P2 J \p\ 2 V V\V2 P 2 J 

p®p 

We choose /3 > and assume for simplicity that the curvature term k is the standard mean 
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curvature : 



k(x) = div 



rv¥l 



d 2 <s> 

ftr 2 dy 2 



|V$| 



' 9x <9j/ dxdy 



X) 

9?/ / dy 2 



V 



Tr(D 2 $(x)) 



|V*r 1 Tr 



|V$| 2 



/V$(x) (8) V$(x) 



(x) 



|V*| 1 ' 1 V |V$(x)| 2 

where D 2 & is the Hessian matrix of <!>. Equation (|3.15p becomes 
^(t,x)+g(|V/(x)| 2 )F(x)|V$(x)| 



£> 2 $(x) 



5 (|V/(x)| 2 )/? Tr(Z) 2 cI>(x))-Tr 



V$(x) (8) V$(x) 



|V$(x) 



D 2 $(x) 



So 



+ 5 (|V/(x)| 2 )(l - |^(x)|)(y(x), V*(x)> = 
$ t (i,x)+ 5 (|V/(x)| 2 )F(x)|V<J?(x)| 



V$(x) V$(xl 



|V$(x) 



-/?g(|V/(x)| 2 )Tr 
+ 5 (|V/(x)| 2 )(l - |#(x)|)(F(x), V*(x)> = 
where I is the identity matrix. Setting 



D 2 $(x 



A(p) := I 



and F(pc,p, X) :- 



(4.19) 



(4.20) 



(4.21) 



(4.22) 



<7(|V/(x)| 2 )tf (x)|p| - /?<?(|V/(x)| 2 ) Tr (A(p)X) + g(\ V/(x)| 2 )(l - |tf (x)|)<F(x),p) (4.23) 
we get 

+ F(x,V$,L> 2 $) = . (4.24) 

The Hamiltonian F : I 2 x I 2 x 52 — ► R, is independent of t and Here 5/v is N x N the 

symmetric matrices space endowed with th classical order. 

First, we have to verify we may use the viscosity theory framework. 
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Definition 4.1 Let T be a function : 



x I 2 x R 2 x ^2 — » M 
(t,a;,p,X) 1 ^ F(t,x,p,X) 



T is proper if : 

- it is a non- decreasing function with respect to the first variable t : 

T(t,x,p,X) < F(a,x,p,X) ift < s. (4.25) 

- it is a decreasing function with respect to the last variable X 

T(t, x, p, X) < T(t, x, p, Y) siY < X. (4.26) 



Theorem 4.1 The function F is proper 

Proof - Since F does not explicitly depend on t, F(t,x,p, X) = F(x,p, X), the condition 
(|4.25p is verified. 

Assuming p / 0. The matrix A{p) defined by (I4.22h is semi-definite positive, and can be 
written as A = a a 1 . So we have 

2 

Tr (AX) = Tr (a^X) = Tr (a l Xa) = ^ a\Xa u 

i=i 

where dj stands for the ith column of a. 

Let X,Y £ S2 be given. Let us assume Y > X, then : 

Vt G {1,2}, alXo-i < a\Yo-i. 

The constant (3 is positive and the function g as well, so we deduce that : 

-/3 9 (|V/(x)| 2 )Tr(^(p)X) > -f3g(\VI(x)\ 2 )Tr(A(p)Y). 

so F(x,p,X)<F(x t p,Y). 

We conclude that the function F is proper. □ 
The general viscosity theory framework is thus well posed. 

4.2 Ishii and Sato theorem 

In what follows we use a theorem by Ishii and Sato [6] that gives existence of viscosity 
solution to singular degenerate parabolic (Hamilton-Jacobi) equations with nonlinear oblique 
derivative boundary conditions. 

In this subsection we recall this theorem. In the sequel we denote p the function defined 

from E 2 in R by p(p, q) = min ( — ^ ^) 1 J . 

V mm (M>M) / 
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Theorem 4.2 ([6] Theorem 2.1 p 1079 ) 
HI- f GC([0,T] xTixRx (R 2 \{0 r , 2 }) xS 2 ), 

H2- There exists 7 G R suc/i that for every (t,x,p,X) e[0,T]xOx (R 2 \ {0r 2 }) x S 2 , the 
function A 1— > x, A,p, X) — 7A is non- decreasing on R. 

H3— For every R > 0, there exists a continuous non- decreasing function ipn : [0, +oo[— ► 
[0, +oo[such that </3r(0) = and for every X,Y £ S 2 and fii,fi 2 6 [0, +oo[ satisfying : 

I 0\ ( I —l\ (I 

x y J^V/ / ) +fl2 {o 1 

then 

F(t, x, A, p, X) - F(t, y, A, q, -Y) > 

/ 9 9 \ \ ' / 

-¥>fl(MHF-y| + / 9 (P>#) ) + |p- g| + \x - y\(l + max(|p|, \q\))) , 

Vt G [0,T], x,y G H, A G R suc/i f/iat |A| < i? and p, g G R 2 \ {0 R2 }. 

iJ4— J 7 is continuous at (t, x, A, 0, 0) pour tout (t, x, A) G [0, T] x $7 x R in the following sense: 

-00 < F*(t,x,\,0,Q) = J*(t,x,A,0,0) < +00 

where T* (respectively T* ) are the upper (respectively lower) semi- continuous envelopes 
of T, defined on [0, T] x Q x R x R 2 x S 2 . 

Bl- B G C(R 2 x R 2 ) n C 1 - 1 ^ 2 x (R 2 \ {0 R2 })). 

B2— Pour tout x G R 2 , the function p ^ B(x,p) 1-positively homogeneous with respect to p, 
i.e., B(x, Xp) = \B(x,p),\/\ > 0,p G R 2 \ {0 R2 }. 

B3— There exists a positive constant 9 such that (v(z), D p B(z,p)) > 9 for every z G d£l and 
pGR 2 \{0 R2 }. 

Here v{z) is the unit outer normal vector of Q at z G dQ. 

Assume [HI, H2, H3, HA, Bl, B2, B3] are satisfied and consider the following problem 

/ cp t + H^x )( p,V<t>,D 2 cp) =0, m Q :=]0,T[xn 
{) {B(x,V(l)) = on]0,T[xdn { > 

• Let (f> G USC([0, T[x£l) and ip G LSC([0,T[xQ) be, respectively, viscosity sub and super- 
solutions of U.28\) . If (f)(0, x) < ip(0, x) for x G H, then <j)<ip sur ]0, T[xU. 

• For every function g G C(0) there exists a unique viscosity solution <j) G C([0,T[xQ) of 
{4. 28 ) such that 4>(0,x) = g{x) on £1. 

Here USC([0, T[xQ) (respectively LSC([0, T[xfi) denote the set of upper semicontinuous 
(respectively lower semicontinuous) functions. 
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4.3 Existence and uniqueness of the solution of the evolution problem 

In the sequel we assume that the balloon force H verifies (H2) '■ 

(TC2) '■ H is Lipschitz continuous on S7. 

In order to use Ishii-Sato theorem we have to verify every hypothesis. The Hamiltonian 
F(t,x,r,p,X) = F(x,p,X) is defined by (1031 : 

F(x,p,X) = £(x)tf(x)bp| -Pg(x)Tr(A(p)X) +flf(x)(l - \H(x)\)(V(x),p) . (4.29) 
Let us define the symmetric, semi-definite positive matrix A(x,p) : 

A( X ,p)=Pg(x)(I-^f), 
\P\ 

so that 

F(x, p,X) = - Tr (A(x, p)X) + g(x)H(x) \p\ + £(x) (1 - \H(y:)\)(V(x.),p) . (4.30) 

In this case F does not depend neither on t nor A. 
We choose a Neumann-type boundary condition : 

(v^(x), i/(x)) K2 = |^(x) = on r = an 

[HI] Function F G C(f2 x (R 2 \ {Or 2 }) x S 2 - p = is a singular point. 

Theorem 12.31 yields that the gradient vector flow V is continuous on [0, T] x Q. 

In addition we assumed that the balloon function H is continuous on Q . 
Therefore, the Hamiltonian F is continuous on x R 2 \ {0^2} x £2 

[H2] Let us show there exists a constant 7 6 R such that for each (x,p, X) G xfi x (R 2 \ 
{0 K 2}) x 52, the function A 1— > F(x,p,X) — 7A is non-decreasing on R. Since F does 
not explicitly depend on A, any negative constant 7 is suitable. 

[H3] As F does not depend on t and r, we have to find a continuous increasing function 
cp : [0,+oo[— > [0, +oo[ satisfying ip(0) = such that if X, Y G 1S2 and /ii,//2 £ [0, +oo[ 
satisfy : 



then Vx, y G fi, and p,ggK 2 \ 

F(x,p,X)-F(y, g ,-Y) 

> -99 (/ii(|x - y| 2 +p(p,q) 2 ) +fi 2 + \p~q\ + |x- y|(l +max( 
Let use the following lemma [8] : 



4 PROPAGATION EQUATION STUDY 



16 



Lemma 4.1 Ifp,q£R N \ {0 R2 }, then : 



_p g_ 

W \q\ 



^ P-Q i \ 

mm p , a 



Given X, 1" G £2 an d Mis S [0, +00 [ verifying (|4.31|) . Let be r, s G K 2 , so we have 
(Xr,r) + (Ys,s) < m\r - s\ 2 + fi 2 (\r\ 2 + |s| 2 ) 

Let x, y G $7 and p,g 6 K 2 \ {0r2}. Following C. Le Guyader [8], we split F in three 
terms and verify [H3] for each term. 
F(x.,p,X)-F(y,q,-Y) = 

-(Tr(A(x,p)X)+Tr(A(y,g)r)) 

V v ' 

(a) 

+ g(x)ff(x)|p|-g(y)ff(y)|g| 

(6) 

+ 5(x)(l- |H(x)|)(V(x),p) M2 -s(y)(l- | J e"(y)|)(T>(y),g) R2 . 

V v ' 

(c) 

• (a) estimate.- As A(x,p) = <j(x,p)cr*(x,p), 

Tr(A(x,p)X)+Tr(,4(y,g)Y) < 

Ml Tr ( (<r(x, p) - a(y, q) ) (<r(x, p) - a(y, q) )*) (Tr (<r(x, (y, p)) + Tr ((r(y, g)^ (y, q))) 
< MiTr ((<r(x,p) - <r(y,g))(<r(x,p) - <r(y,g))') +2^ 2 /35 

V V ' 

(al) 

where <5 is an upper bound of 5 on O (for example ). Expression (al) verifies : 

v27rcr 



Tr ((cr(x,p) - a(y,q))(a(x,p) - a{y,q)Y) 

= Tr (o-(x,p)<7*(x,p) - o-(x,p)<T t (y,g) - cr(y, g)cr*(x,p) + cr(y, g)cr* (y, g)) 
= Tr {A(x,p) -fj(x,p)(T t (y,g) - <r(y, g)a'(x,p) + A(y,q)) 

= 0g(x) + 2fiyJW)^W) (-2 + |^j(Pi5i +^g 2 )l + /9£(y) 
and with 

-2 + ^|(Pi«i+P2<6)) <-j^(Pi«i +!»*.) 
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wc get : 



Tr ((ff(x,p) - ff(y,g))((j(x,p) - <j(y,g)) f ) 



< 0£(x) - 2/3y^(x) y^(y)_(p igi + p 2 g 2 ) + #?(y) 



We finally obtain 

Tr (A(x, p)X) + Tr (A(y, g)Y) < /ii 



Moreover we have the relation : 



So we get by lemma 13.11 : 



< 



( - Vffly)) I, + v^y) ( 



+ 2/x 2 /M (4.33) 



P 



<2/3(yjWj 



\\p\ \q\ 

2 

P Q 



<2/3iY 1 2 |x-y| 2 + 2/35 



P Q_ 

p\ \q\ 



p_ q_ 

\P\ \q\ 



Thus we conclude with lemma [4TT1 



Tr (A(x,p)X) + Tr (A(y, q)Y) < Ml (2/3K 2 |x - y| 2 + 2/W4p(p, q) 2 ) + 2^(35 (4.34) 

• (b) estimate.- We have assumed H to be Lipschitz continuous on 0, so x 6 fi i— > 
g(x)H(x) is Lipschitz continuous and bounded as well. Let be x, y G f2 

|flf(x)ff(x)|p| -g(y)H(y)\ q \\ < \(g(x)H(x)-g(y)H(y)) \p\\ + \g(y)H(y)\\p\ - \q\\\ 

< K 2 |x-y|max(|p|,|g|)+0||p| - \q\\ 

< K 2 \x - y\ max(|p|, \q\) + 9\p - q\ 

where K 2 is the Lipschitz-constant of the function x £ $7 i— > g(x)H(x) and 9 is a bound 

ofxeH^ |g(x)ff(x)|. 
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• (c) estimate.- We have assumed (TC2) sox£(1h g(x)(l — |iJ(x)|)1^(x) is Lipschitz 
continuous and bounded. 



Sf(x)(l - |F(x)|)(y(x),p) R2 -g(y)(l - \H(y)\)(V(y),q) 



< 



<5(x)(l - |i?(x)|)y(x) -g(y)(l - \H( y )\)V(y),ph 



+ 



(g(y)(l-\H(y)\)V(y),p-q) 1 



< K 3 \x - y| max(|p|, \q\) + (\p - q\ 

where K 3 is the Lipschitz-constant of the function x G ^ g(-x)(l — \H(x)\)V(x) and 
C is a bound of x e H ^ |flf(x)(l - |F(x)|)t>(x)|. 
Finally relation (|4.32j) gives 

- (F(x,p, X) - F(y, g, -Y)) < fit (2/3# 2 |x - y| 2 + 8/Mp(p, g) 2 ) + 2p 2 (35 

+ (K 2 |x - y| max(|p|, \q\) + 6>|p - g|) 
+ (K 3 |x - y| max(|p|, |g|) + (\p - q\). 

< max(2/^ 2 , 8/M, K 2 +K 3 ,9 + C)[(MI* - y| 2 + pip, q) 2 ) 
+fi 2 ) + max(\p\, |<?|)|x - y| + |p - </|] 

< max(2/3^ 2 , 8(35, K 2 + K 3 ,6 + Q ( W (|x - y| 2 + p(p, q) 2 ) 
+ ii2 + (1 + max(|p|, M))|x - y| + \p- q\) . 

Finally 

F{-K,p,X)-F{y,q,-Y) 

> -ip (mi(|x - y| 2 + p(p, q) 2 ) + ll 2 + (1 + max(|p|, |g|))|x - y| + \p - q\)) 
where the function ip is defined by : 

ip(m) = max(2/3Kf, 8(35, K 2 + K 3 ,6 + Qm 
Hypothesis [H3] is then verified. 
[H4] F is continuous at (x, 0, 0) for any x£!l because F*(x, 0,0) = F*(x, 0, 0) = 0. 



For the three last hypothesis Bl, B2 and B3 on the boundary condition, the proof is 
the same as in C. Le Guyader [8]. 
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[Bl] This hypothesis consists in showing that the function B which defines the boundary 
condition on ]0, +oo[x<9ft is C(R N x R N ) n C 1 - 1 ^ x (R N \ {0 RN })). We have chosen 
a Neumann-type condition, by denoting v (x) the outward unit normal vector of dO, at 
the point x, our boundary condition is written : 

£(x,p) = (z/(x),p) R2 . 

The point [SI] will be satisfied if v is a C ' vector field which is the case if f2 has a C 2 
boundary. 

[B2] B is 1-positively homogeneous with respect to p : 

B(x, Xp) = (u(x), Ap) R2 = AB(x,p), VA > 0,p G M 2 \ {0 R2 }. 

[B3] Let zedn. 

(u(z),D p B(z,p)) R 2 = \u(z)\ 2 = 1. 
The last condition is verified with 6 = 1, 

We may now conclude since assumptions [ HI, H2, H3, HA, Bl, B2 , B3] are satisfied. 

Theorem 4.3 Assume that the image function I verifies the hypothesis (TCi) and the ballon 
force H satisfies (Ti.2)- Consider the following problem 



* t (t,x) - 5 (|V/(x)| 2 ) ((/3 K (x) - ff(x))|V*(x)| - (1 - |tf(x)|)(F(x), V$(x))J = 

in ]0, +oo[xf2, 

— (x) = on \0,+oo[xdQ , 

(4-35) 

• Let $ 6 {/SCQO, T[x f2) and ^ G LSC([0,T[xQ) be, respectively, viscosity sub and super- 
solutions of : If®(Q,x) < ~®(Q,x) for x G IT, t/ien $ < ^ in ]0,T[xIT. 

• For every g G C(fi), i/tere is a unique viscosity solution <I> G C([0,T[xfi) o/ satisfying 
3>(0,x) = ^(x) on $7. 

• Equation k3.15\) has a unique viscosity solution $ G C([0,T[xO). 



5 Conclusion 

We have recalled Gradient Vector Flow model that we hahe justified and we have given 
regularity properties. Then we proved existence and uniqueness fo viscosity solutions of the 
Hamilton-Jacobi equation derived fron the GVF-geodesic active contour model. 

Next step is to perform the numerical realization of this GVF-geodesic active contour 
process solving he Hamilton-Jacobi equation (I4.35p . We shall use his method to the perform 
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"tuffeau" tomographic images segmentation. This material is has been used during past 
centuries to build monuments as castles and churches in the Val-de-Loire area. These images 
allow to get information on the structure of the damaged material : we have to identify 
different phases as calcite (light grey), silice (dark grey) and porosity (black). 

We shall combine the GVF-geodesic model with a region segmentaion approach to identify 
the three constituents of tuffeau . The segmentation of these images is a step of pretreatment 
which aims at reconstructing the stone porosity domain. 
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